function [model_settings,info] = model_id(model_settings0,w,data)
%MODEL_ID identify nu and sigma2 according to meanR2 or varR2 or meanM2 using
%  nonlinear least square fit.
%
% [model_settings] = model_id(model_settings0) returns the default models as Gangshi
% identified.
%
% [model_settings] = model_id(model_settings0,w,data) identified using data
%   specified by DATA.
%
% INPUTS:
% w: a array
%
% OUTPUTS:
% model_settings = a structure get from MODEL_ID().
% 
% See also: model.m

model_settings = model_settings0;

model_settings0.W       = [0.10,0.15,0.20,0.50,1.0];
model_settings0.nu      = [0.0451e-4,1.3915e-05,0.3583e-4,0.1232e-3,0.3939e-3];
model_settings0.sigma2  = [1.0810e-3,6.2341e-3,2.1983e-2,0.0990,0.3007];
model_settings0.Rh      = [0.1002,0.1509,0.2021,0.5129,1.0432];
model_settings0.K       = [0.99962,0.9955,0.9914,0.9730,0.9532];
model_settings0.Tau     = [10.0740,7.6978,5.3215,2.1915,1.1525];

% If no data is provided, use the default model_settings
if nargin == 1
    model_settings.W      = model_settings0.W;
    model_settings.nu     = model_settings0.nu;
    model_settings.sigma2 = model_settings0.sigma2;
    model_settings.Rh     = model_settings0.Rh;
    model_settings.K      = model_settings0.K;
    model_settings.Tau    = model_settings0.Tau;
else
    model_settings.W = w;
    nlinfitoptions = statset('Display','off');
    for i = 1:length(w)
        fprintf(1,'Processing #%d data: w = %f\n',i,w(i));
        t = data(i).t;
        indx = (t<=model_settings0.FittingTimeRange);

        % Fitting of Rh, K and Tau
        % FIXME: Need to update this part.
        if isempty(data(i).h)
            Rh  = interp1(model_settings0.W,model_settings0.Rh,w(i));
            K   = interp1(model_settings0.W,model_settings0.K,w(i));
            Tau = interp1(model_settings0.W,model_settings0.Tau,w(i));
        end
        model_settings.K(i)   = K;
        model_settings.Rh(i)  = Rh;
        model_settings.Tau(i) = Tau;

        % Fitting of nu and sigma2
        beta0(1) = interp1(model_settings0.W,model_settings0.nu,w(i));
        beta0(2) = interp1(model_settings0.W,model_settings0.sigma2,w(i));
        switch(model_settings0.FitTo)
            case 'meanR2'
                meanR2 = data(i).meanR2;
                [beta] = nlinfit(t(indx),meanR2(indx),@model_meanR2,beta0,nlinfitoptions);
                %             meanR2_pred(:,i) = model_meanR2(beta,t);
            case 'VarR2'
                varR2 = data(i).varR2;
                [beta] = nlinfit(t(indx),varR2(indx),@model_varR2,beta0,nlinfitoptions);
            case 'meanM2'
                meanSlope = data(i).meanM2;
                [beta] = nlinfit(t(indx),meanSlope(indx),@model_meanSlope,beta0,nlinfitoptions);
            otherwise
                error('LPCVD:model_id:Bad model.FitTo.');
        end
        model_settings.nu(i)  = beta(1);
        model_settings.sigma2(i) = beta(2);
    end

    % Plot the identified models and compare with data
    for i = 1:length(w)
        t  = data(i).t;
        dt = t(2)-t(1);
        h       = zeros(length(t),1);
        rho     = zeros(length(t),1);
        meanR2  = zeros(length(t),1);
        varR2   = zeros(length(t),1);
        State.h   = 0;
        State.rho = 0;
        State.meanAlpha2 = zeros(1,20);
        State.meanBeta2  = zeros(1,20);
        State.meanR2     = 0;
        State.varAlpha2  = zeros(1,20);
        State.varBeta2   = zeros(1,20);
        State.varR2      = 0;
        for j = 2:length(t)
            State     = model(w(i),State,dt,model_settings);
            h(j)      = State.h;
            rho(j)    = State.rho;
            meanR2(j) = State.meanR2;
            varR2(j)  = State.varR2;
        end
        msg_title = sprintf('w = %f',model_settings.W(i));

        figure;
        subplot(2,1,1);
        meanR2_pred = model_meanR2([model_settings.nu(i),model_settings.sigma2(i)],t);
        plot(t,data(i).meanR2,t,meanR2,t,meanR2_pred);
        xlim([0,options.t]);
        ylabel('mean(R2)');
        title(msg_title);

        subplot(2,1,2);
        varR2_pred = model_varR2([model_settings.nu(i),model_settings.sigma2(i)],t);
        plot(t,data(i).varR2,t,varR2,t,varR2_pred);
        ylabel('var(R2)');
        xlim([0,options.t]);
    end
end


end
%%
function meanR2_model = model_meanR2(beta,t)
nu = beta(1);
sigma2 = beta(2);

meanR2_model = zeros(length(t),1);
for i = 1:length(t)
    meanAlpha2 = zeros(1,20);
    meanBeta2  = zeros(1,20);
    for j = 1:20
        temp  = sigma2/(2*nu*j^2);
        temp2 = exp(-2*nu*t(i)*j^2);
        meanAlpha2(j) = temp+(meanAlpha2(j)-temp)*temp2;
        meanBeta2(j)  = temp+(meanBeta2(j)-temp)*temp2;
    end
    meanR2_model(i) = sum(meanAlpha2+meanBeta2)/(2*pi);
end
end

%%
function varR2 = model_varR2(beta,t)
nu = beta(1);
sigma2 = beta(2);

varR2 = zeros(length(t),1);
for i = 1:length(t)
    meanAlpha2 = zeros(1,20);
    meanBeta2  = zeros(1,20);
    varAlpha2  = zeros(1,20);
    varBeta2   = zeros(1,20);
    for j = 1:20
        temp1 = exp(-2*nu*j^2*t(i));
        temp2 = exp(-4*nu*j^2*t(i));
        temp3 = sigma2*(temp1-1)/(-2*nu*j^2);
        varAlpha2(j) = temp2*varAlpha2(j)+4*temp1*temp3*meanAlpha2(j)+2*temp3^2;
        varBeta2(j)  = temp2*varBeta2(j)+4*temp1*temp3*meanBeta2(j)+2*temp3^2;
    end
    varR2(i) = sum(varAlpha2+varBeta2)/(4*pi^2);
end
end

%%
function meanM2 = model_meanSlope(beta,t)
nu = beta(1);
sigma2 = beta(2);

meanM2 = zeros(length(t),1);
for i = 1:length(t)
    meanAlpha2 = zeros(1,model_settings0.mode);
    meanBeta2  = zeros(1,model_settings0.mode);
    for j = 1:model_settings.mode
        temp  = sigma2/(2*nu*j^2);
        temp2 = exp(-2*nu*t(i)*j^2);
        meanAlpha2(j) = temp+(meanAlpha2(j)-temp)*temp2;
        meanBeta2(j)  = temp+(meanBeta2(j)-temp)*temp2;
    end
    meanM2(i) = sum(meanAlpha2.*model_settings0.M2ModeWeighting(:,1)+meanBeta2.*model_settings0.M2ModeWeighting(:,2));
end
end

